Models for Q2

Erick Calderon-Morales and Leland Werden


library(dplyr)
library(nlme)
library(purrr)
library(performance)
library(tibble)
library(here)
library(reactablefmtr)
library(reactable)
library(emmeans)
library(car)
library(ggplot2)
library(modelr)
library(knitr)

Load functions and data

# Load all joined dataset
source("./scripts/code_join_data_full_dataset.R")

# Step was done like this because I am working with a subset of the data
# source cleaned data
source("./scripts/code_clean_data_nodules.R")
# Load custom made functions
source("./R/functions_models.R")
source("./R/function_plots.R")
source("./R/function_validation_plots.R")
source("./R/function_for_inference_emmeans_and_percentage_diff.R")

Q2: How does increased nutrient and/or water availability influence seedling water- and nutrient-use traits and the relationships with N-fixing bacteria?

Model fitting

Traits

# Take response variables names
response_vars_q2 <-
  set_names(c("amax", "gs", "sla"))
## Create empty list
models_q2 <- list()
model_q2_amax <- lme(amax ~ nfixer*treatment +init_height,
                                        random = ~1|spcode,
                                        data = data_for_models)

model_q2_amax <- list(model_q2_amax)

names(model_q2_amax) <- "amax"
model_q2_gs <- lme(gs ~ nfixer*treatment +init_height,
                                        random = ~1|spcode,
                                        data = data_for_models)

model_q2_gs <- list(model_q2_gs)

names(model_q2_gs) <- "gs"
model_q2_sla <- lme(sla ~ nfixer*treatment +init_height,
                                        random = ~1|spcode,
                                        data = data_for_models)

model_q2_sla <- list(model_q2_sla)

names(model_q2_sla) <- "sla"                      
## PNUE
model_q2_pnue_log <- lme(log(pnue) ~ nfixer*treatment +init_height,
                                        random = ~1|spcode,
                                        data = data_for_models)

model_q2_pnue_log <- list(model_q2_pnue_log)

names(model_q2_pnue_log) <- "pnue_log"
model_q2_pnue_nlme <- lme(pnue ~ nfixer*treatment + init_height,
                                        random = ~1|spcode,
                                        weights = varIdent(form = ~1|spcode*treatment),
                                        data = data_for_models)

model_q2_pnue_nlme <- list(model_q2_pnue_nlme)

names(model_q2_pnue_nlme) <- "pnue_nlme"
## Narea_g_m2 log model
model_q2_n_area_log <- lme(log(narea_g_m2) ~ nfixer*treatment + init_height,
                                                random = ~1|spcode,
                                                data = data_for_models)

model_q2_n_area_log <- list(model_q2_n_area_log)

names(model_q2_n_area_log) <- "n_area_log"
model_q2_n_area_nlme <- lme(narea_g_m2 ~ nfixer*treatment + init_height,
                                        random = ~1|spcode,
                                        weights = varIdent(form = ~1|spcode*treatment),
                                        data = data_for_models)

model_q2_n_area_nlme <- list(model_q2_n_area_nlme)

names(model_q2_n_area_nlme) <- "n_area_nlme"
## WUE log model
model_q2_wue_log <- lme(log(wue) ~ nfixer*treatment + init_height,
                                        random = ~1|spcode,
                                        data = data_for_models)

model_q2_wue_log <- list(model_q2_wue_log)

names(model_q2_wue_log) <- "wue_log"
model_q2_wue_nlme <- lme(wue ~ nfixer*treatment + init_height,
                                        random = ~1|spcode,
                                        weights = varIdent(form = ~1|spcode*treatment),
                                        data = data_for_models)

model_q2_wue_nlme <- list(model_q2_wue_nlme)

names(model_q2_wue_nlme) <- "wue_nlme"

Nodule colonization

# Delete unused variables
data_nodules_cleaned <-
    data_nodules_cleaned %>%

        # add id to rownames for keep track of the rows
        column_to_rownames("id") %>%
        dplyr::select(spcode, treatment, everything())

Nodule weight

nlme_nodule_weight <- lme(estimated_total_nodule_mass_per_plant ~ treatment + init_height,
                                    random = ~1|spcode,
                                    weights = varIdent(form = ~1|spcode*treatment),
                                    data = data_nodules_cleaned)


model_q2_nodule_weight <- list(nlme_nodule_weight)

names(model_q2_nodule_weight) <- "nodule_weight"

Nodule count

nlme_nodule_count <- lme(total_number_of_plant_nodules ~ treatment + init_height,
                                    random = ~1|spcode,
                                    weights = varIdent(form = ~1|spcode*treatment),
                                    data = data_nodules_cleaned)


model_q2_nodule_count <- list(nlme_nodule_count)

names(model_q2_nodule_count) <- "nodule_count"
# Append models to model list
models_q2 <- append(model_q2_amax, models_q2)
models_q2 <- append(model_q2_gs, models_q2)
models_q2 <- append(model_q2_sla, models_q2)

models_q2 <- append(model_q2_n_area_log, models_q2)
models_q2 <- append(model_q2_n_area_nlme, models_q2)

models_q2 <- append(model_q2_pnue_log, models_q2)
models_q2 <- append(model_q2_pnue_nlme, models_q2)

models_q2 <- append(model_q2_wue_log, models_q2)
models_q2 <- append(model_q2_wue_nlme, models_q2)

models_q2 <- append(model_q2_nodule_count, models_q2)
models_q2 <- append(model_q2_nodule_weight, models_q2)
names(models_q2)
##  [1] "nodule_weight" "nodule_count"  "wue_nlme"      "wue_log"      
##  [5] "pnue_nlme"     "pnue_log"      "n_area_nlme"   "n_area_log"   
##  [9] "sla"           "gs"            "amax"

Model Assumptions

Maximal photosynthesis

validation_plots(models_q2$amax, data = data_for_models, group = "spcode")

Stomatal Conductance

validation_plots(models_q2$gs, data = data_for_models, group = "spcode")

SLA

validation_plots(models_q2$sla, data = data_for_models, group = "spcode")

Water Use Efficiency

validation_plots(models_q2$wue_log, data = data_for_models, group = "spcode")

validation_plots(models_q2$wue_nlme, data = data_for_models, group = "spcode")

PNUE

validation_plots(models_q2$pnue_log, data = data_for_models, group = "spcode")

validation_plots(models_q2$pnue_nlme, data = data_for_models, group = "spcode")

Nitrogen concentration per unit of area

validation_plots(models_q2$n_area_log, data = data_for_models, group = "spcode")

validation_plots(models_q2$n_area_nlme, data = data_for_models, group = "spcode")

Nodule weight

validation_plots(models_q2$nodule_weight, data = data_nodules_cleaned, group = "treatment")

Nodule count

validation_plots(models_q2$nodule_count, data = data_nodules_cleaned, group = 'spcode')

Model inference

## r2 models
models_q2 %>%
    map(., r2) %>%
    unlist()
## Random effect variances not available. Returned R2 does not account for random effects.
## nodule_weight.R2_conditional.Conditional R2 
##                                     0.41047 
##       nodule_weight.R2_marginal.Marginal R2 
##                                     0.13936 
##  nodule_count.R2_conditional.Conditional R2 
##                                     0.95127 
##        nodule_count.R2_marginal.Marginal R2 
##                                     0.04775 
##      wue_nlme.R2_conditional.Conditional R2 
##                                     0.90104 
##            wue_nlme.R2_marginal.Marginal R2 
##                                     0.83953 
##       wue_log.R2_conditional.Conditional R2 
##                                     0.71998 
##             wue_log.R2_marginal.Marginal R2 
##                                     0.55827 
##     pnue_nlme.R2_conditional.Conditional R2 
##                                     0.54153 
##           pnue_nlme.R2_marginal.Marginal R2 
##                                     0.25717 
##      pnue_log.R2_conditional.Conditional R2 
##                                     0.49405 
##            pnue_log.R2_marginal.Marginal R2 
##                                     0.19552 
##   n_area_nlme.R2_conditional.Conditional R2 
##                                     0.85290 
##         n_area_nlme.R2_marginal.Marginal R2 
##                                     0.50938 
##    n_area_log.R2_conditional.Conditional R2 
##                                     0.86656 
##          n_area_log.R2_marginal.Marginal R2 
##                                     0.58168 
##           sla.R2_conditional.Conditional R2 
##                                          NA 
##                 sla.R2_marginal.Marginal R2 
##                                     0.28831 
##            gs.R2_conditional.Conditional R2 
##                                     0.46869 
##                  gs.R2_marginal.Marginal R2 
##                                     0.39313 
##          amax.R2_conditional.Conditional R2 
##                                     0.94110 
##                amax.R2_marginal.Marginal R2 
##                                     0.49635
## r2 models
models_q2 %>%
    map(., AIC) %>%
    unlist()
## nodule_weight  nodule_count      wue_nlme       wue_log     pnue_nlme 
##       -65.217       411.462       302.027       124.712      1048.375 
##      pnue_log   n_area_nlme    n_area_log           sla            gs 
##        85.397        88.499        -2.973      -925.487      -159.109 
##          amax 
##       483.009

Anova tables

#map(models_q2, ~Anova(.x, type = "III", test.statistic = c("F")))
map(models_q2, ~anova.lme(.x, type = "marginal"))
## $nodule_weight
##             numDF denDF F-value p-value
## (Intercept)     1    41   3.298  0.0767
## treatment       3    41   4.623  0.0071
## init_height     1    41   1.376  0.2476
## 
## $nodule_count
##             numDF denDF F-value p-value
## (Intercept)     1    41   3.017  0.0899
## treatment       3    41  13.952  <.0001
## init_height     1    41   0.732  0.3971
## 
## $wue_nlme
##                  numDF denDF F-value p-value
## (Intercept)          1   102  121.59  <.0001
## nfixer               1     6    0.34  0.5798
## treatment            3   102   66.94  <.0001
## init_height          1   102    0.63  0.4307
## nfixer:treatment     3   102   49.70  <.0001
## 
## $wue_log
##                  numDF denDF F-value p-value
## (Intercept)          1   102   36.36  <.0001
## nfixer               1     6    4.09  0.0897
## treatment            3   102   26.88  <.0001
## init_height          1   102    0.40  0.5264
## nfixer:treatment     3   102    4.67  0.0042
## 
## $pnue_nlme
##                  numDF denDF F-value p-value
## (Intercept)          1   102   89.71  <.0001
## nfixer               1     6    0.41  0.5455
## treatment            3   102   16.00  <.0001
## init_height          1   102    7.83  0.0061
## nfixer:treatment     3   102    6.70  0.0004
## 
## $pnue_log
##                  numDF denDF F-value p-value
## (Intercept)          1   102   947.5  <.0001
## nfixer               1     6     0.5  0.4992
## treatment            3   102     3.6  0.0156
## init_height          1   102     2.2  0.1425
## nfixer:treatment     3   102     1.8  0.1541
## 
## $n_area_nlme
##                  numDF denDF F-value p-value
## (Intercept)          1   102  15.472  0.0002
## nfixer               1     6  11.956  0.0135
## treatment            3   102  20.985  <.0001
## init_height          1   102   1.109  0.2947
## nfixer:treatment     3   102   2.026  0.1149
## 
## $n_area_log
##                  numDF denDF F-value p-value
## (Intercept)          1   102   0.606  0.4381
## nfixer               1     6  15.964  0.0072
## treatment            3   102   7.372  0.0002
## init_height          1   102   0.575  0.4501
## nfixer:treatment     3   102   2.385  0.0735
## 
## $sla
##                  numDF denDF F-value p-value
## (Intercept)          1   102  118.25  <.0001
## nfixer               1     6    0.43  0.5342
## treatment            3   102    0.25  0.8618
## init_height          1   102    6.32  0.0135
## nfixer:treatment     3   102    1.80  0.1521
## 
## $gs
##                  numDF denDF F-value p-value
## (Intercept)          1   102   8.874  0.0036
## nfixer               1     6   6.401  0.0447
## treatment            3   102  18.662  <.0001
## init_height          1   102   1.015  0.3161
## nfixer:treatment     3   102   1.698  0.1722
## 
## $amax
##                  numDF denDF F-value p-value
## (Intercept)          1   102  12.364  0.0007
## nfixer               1     6   7.969  0.0302
## treatment            3   102   0.257  0.8562
## init_height          1   102   5.555  0.0203
## nfixer:treatment     3   102   7.384  0.0002

Post-Hoc: Tukey’s test

Maximal photosynthesis

as_tibble(emmeans(models_q2$amax,
        pairwise ~ treatment*nfixer,
        adjust = "tukey"
        )$contrast) %>%
        mutate(across(2:6, round, 6)) %>%
        kable()
contrast estimate SE df t.ratio p.value
no_additions nonfixer - plus_nutrients nonfixer -0.1287 0.5156 102 -0.2495 1.0000
no_additions nonfixer - plus_water nonfixer -0.4310 0.5208 102 -0.8277 0.9911
no_additions nonfixer - plus_water_nutrients nonfixer -0.0716 0.4988 102 -0.1436 1.0000
no_additions nonfixer - no_additions fixer -9.0888 3.2195 6 -2.8230 0.2419
no_additions nonfixer - plus_nutrients fixer -7.0199 3.2267 6 -2.1756 0.4646
no_additions nonfixer - plus_water fixer -11.4230 3.2202 6 -3.5473 0.1120
no_additions nonfixer - plus_water_nutrients fixer -9.7834 3.2227 6 -3.0357 0.1931
plus_nutrients nonfixer - plus_water nonfixer -0.3024 0.5379 102 -0.5621 0.9992
plus_nutrients nonfixer - plus_water_nutrients nonfixer 0.0570 0.5175 102 0.1102 1.0000
plus_nutrients nonfixer - no_additions fixer -8.9601 3.2232 6 -2.7799 0.2532
plus_nutrients nonfixer - plus_nutrients fixer -6.8913 3.2296 6 -2.1338 0.4829
plus_nutrients nonfixer - plus_water fixer -11.2943 3.2243 6 -3.5029 0.1174
plus_nutrients nonfixer - plus_water_nutrients fixer -9.6548 3.2280 6 -2.9910 0.2025
plus_water nonfixer - plus_water_nutrients nonfixer 0.3594 0.5209 102 0.6900 0.9971
plus_water nonfixer - no_additions fixer -8.6577 3.2230 6 -2.6863 0.2792
plus_water nonfixer - plus_nutrients fixer -6.5889 3.2302 6 -2.0398 0.5257
plus_water nonfixer - plus_water fixer -10.9919 3.2235 6 -3.4099 0.1296
plus_water nonfixer - plus_water_nutrients fixer -9.3524 3.2260 6 -2.8991 0.2232
plus_water_nutrients nonfixer - no_additions fixer -9.0171 3.2193 6 -2.8009 0.2476
plus_water_nutrients nonfixer - plus_nutrients fixer -6.9483 3.2267 6 -2.1534 0.4743
plus_water_nutrients nonfixer - plus_water fixer -11.3513 3.2198 6 -3.5255 0.1146
plus_water_nutrients nonfixer - plus_water_nutrients fixer -9.7118 3.2220 6 -3.0142 0.1975
no_additions fixer - plus_nutrients fixer 2.0688 0.7107 102 2.9108 0.0810
no_additions fixer - plus_water fixer -2.3342 0.6752 102 -3.4573 0.0175
no_additions fixer - plus_water_nutrients fixer -0.6947 0.6816 102 -1.0191 0.9705
plus_nutrients fixer - plus_water fixer -4.4030 0.7131 102 -6.1741 0.0000
plus_nutrients fixer - plus_water_nutrients fixer -2.7635 0.7236 102 -3.8192 0.0055
plus_water fixer - plus_water_nutrients fixer 1.6395 0.6758 102 2.4259 0.2403
# Treatment effects
emmeans_table_tidy(models_q2$amax,
                        formula = "treatment|nfixer",
                        grouping_var = "nfixer")
## [1] "Formula for pairwise comparisons: "
## pairwise ~ treatment | nfixer
## <environment: 0x630b9b0effb0>

Water Use Efficiency

as_tibble(emmeans(models_q2$wue_log,
        pairwise ~ treatment*nfixer,
        adjust = "tukey"
        )$contrast) %>%
        mutate(across(2:6, round, 6)) %>%
        kable() 
contrast estimate SE df t.ratio p.value
no_additions nonfixer - plus_nutrients nonfixer 0.0227 0.1046 102 0.2167 1.0000
no_additions nonfixer - plus_water nonfixer 0.6581 0.1059 102 6.2117 0.0000
no_additions nonfixer - plus_water_nutrients nonfixer 0.6950 0.1015 102 6.8499 0.0000
no_additions nonfixer - no_additions fixer -0.4350 0.2152 6 -2.0215 0.5342
no_additions nonfixer - plus_nutrients fixer -0.6114 0.2197 6 -2.7829 0.2524
no_additions nonfixer - plus_water fixer -0.3098 0.2154 6 -1.4384 0.8145
no_additions nonfixer - plus_water_nutrients fixer -0.2779 0.2164 6 -1.2845 0.8769
plus_nutrients nonfixer - plus_water nonfixer 0.6354 0.1092 102 5.8190 0.0000
plus_nutrients nonfixer - plus_water_nutrients nonfixer 0.6723 0.1049 102 6.4101 0.0000
plus_nutrients nonfixer - no_additions fixer -0.4577 0.2171 6 -2.1083 0.4943
plus_nutrients nonfixer - plus_nutrients fixer -0.6340 0.2213 6 -2.8656 0.2313
plus_nutrients nonfixer - plus_water fixer -0.3325 0.2175 6 -1.5291 0.7735
plus_nutrients nonfixer - plus_water_nutrients fixer -0.3006 0.2189 6 -1.3736 0.8421
plus_water nonfixer - plus_water_nutrients nonfixer 0.0369 0.1059 102 0.3484 1.0000
plus_water nonfixer - no_additions fixer -1.0931 0.2173 6 -5.0308 0.0254
plus_water nonfixer - plus_nutrients fixer -1.2694 0.2218 6 -5.7235 0.0136
plus_water nonfixer - plus_water fixer -0.9679 0.2175 6 -4.4509 0.0443
plus_water nonfixer - plus_water_nutrients fixer -0.9360 0.2183 6 -4.2867 0.0522
plus_water_nutrients nonfixer - no_additions fixer -1.1300 0.2151 6 -5.2535 0.0207
plus_water_nutrients nonfixer - plus_nutrients fixer -1.3063 0.2197 6 -5.9465 0.0113
plus_water_nutrients nonfixer - plus_water fixer -1.0048 0.2153 6 -4.6678 0.0358
plus_water_nutrients nonfixer - plus_water_nutrients fixer -0.9729 0.2161 6 -4.5018 0.0421
no_additions fixer - plus_nutrients fixer -0.1764 0.1445 102 -1.2206 0.9240
no_additions fixer - plus_water fixer 0.1252 0.1372 102 0.9120 0.9843
no_additions fixer - plus_water_nutrients fixer 0.1571 0.1380 102 1.1379 0.9468
plus_nutrients fixer - plus_water fixer 0.3015 0.1448 102 2.0828 0.4330
plus_nutrients fixer - plus_water_nutrients fixer 0.3334 0.1461 102 2.2824 0.3135
plus_water fixer - plus_water_nutrients fixer 0.0319 0.1372 102 0.2325 1.0000
# Treatment effects
emmeans_table_tidy(models_q2$gs,
                        formula = "treatment",
                        )
## [1] "Formula for pairwise comparisons: "
## pairwise ~ treatment
## <environment: 0x630b9afa9080>

PNUE

as_tibble(emmeans(models_q2$pnue_log,
        pairwise ~ treatment,
        adjust = "tukey"
        )$contrast) %>%
        mutate(across(2:6, round, 6)) %>%
        kable()
contrast estimate SE df t.ratio p.value
no_additions - plus_nutrients 0.2100 0.0745 102 2.8193 0.0291
no_additions - plus_water -0.1149 0.0722 102 -1.5899 0.3889
no_additions - plus_water_nutrients 0.0079 0.0715 102 0.1111 0.9995
plus_nutrients - plus_water -0.3248 0.0758 102 -4.2844 0.0002
plus_nutrients - plus_water_nutrients -0.2020 0.0755 102 -2.6757 0.0425
plus_water - plus_water_nutrients 0.1228 0.0722 102 1.6998 0.3291
# Treatment effects
emmeans_table_tidy(models_q2$pnue_log,
                        formula = "treatment",
                        )
## [1] "Formula for pairwise comparisons: "
## pairwise ~ treatment
## <environment: 0x630b98e5f780>

Nitrogen concentration per unit of area

as_tibble(emmeans(models_q2$n_area_log,
        pairwise ~ treatment,
        adjust = "tukey"
        )$contrast) %>%
        mutate(across(2:6, round, 6)) %>%
        kable()
contrast estimate SE df t.ratio p.value
no_additions - plus_nutrients -0.1575 0.0480 102 -3.2823 0.0076
no_additions - plus_water -0.0236 0.0465 102 -0.5074 0.9572
no_additions - plus_water_nutrients -0.0502 0.0461 102 -1.0895 0.6968
plus_nutrients - plus_water 0.1339 0.0489 102 2.7383 0.0361
plus_nutrients - plus_water_nutrients 0.1073 0.0489 102 2.1925 0.1322
plus_water - plus_water_nutrients -0.0266 0.0465 102 -0.5722 0.9402
# Treatment effects
emmeans_table_tidy(models_q2$n_area_log,
                        formula = "treatment",
                        )
## [1] "Formula for pairwise comparisons: "
## pairwise ~ treatment
## <environment: 0x630b9b089510>
as_tibble(emmeans(models_q2$n_area_log,
        pairwise ~ nfixer,
        adjust = "tukey"
        )$contrast) %>%
        mutate(across(2:6, round, 6)) %>%
        kable()
contrast estimate SE df t.ratio p.value
nonfixer - fixer -0.7214 0.1866 6 -3.865 0.0083
# Treatment effects
as.data.frame(emmeans::emmeans(models_q2$n_area_log,
                                specs = pairwise ~nfixer,
                                type = "response",
                                adjust = "tukey")$emmeans) %>%

        janitor::clean_names() %>%
        dplyr::select(response, everything(),
                        # Remove variables
                      -c(df, lower_cl, upper_cl, se)) %>%

        # Rename response to emmean, this is done when models is log
        dplyr::rename_all(funs(stringr::str_replace_all(., "response", "emmean"))) %>%

        # Calculate % difference between control and variable, this assume that
        # that first name is the control

        dplyr::mutate(difference = ((emmean - first(emmean))),
               perc_difference =((emmean - first(emmean) )/first(emmean))*100) %>%

        dplyr::mutate_if(is.numeric, round, 3)
##   emmean   nfixer difference perc_difference
## 1  1.037 nonfixer      0.000             0.0
## 2  2.133    fixer      1.096           105.7

Boxplots traits

# Step done for getting predictions from models for Q2
data_for_predictions <-
    data_for_models %>%

        rownames_to_column("id") %>%

        # Remove unused variables
        dplyr::select(id, spcode, treatment, nfixer, init_height)
string <- c("models_q2")

data_pred_traits <-

        # Get models prediction
        gather_predictions(data_for_predictions ,

                           # Return predictions
                            models_q2$amax,
                            models_q2$sla,
                            models_q2$gs,
                            models_q2$wue_log,
                            models_q2$pnue_log,
                            models_q2$n_area_log

                            ) %>%

        pivot_wider(names_from = model, values_from = pred) %>%
            rename_all(funs(

                # rename columns
                stringr::str_to_lower(.) %>%
                stringr::str_replace(., c(string),"pred_") %>%

                # Remove dollar sing
                gsub("\\$", "", .)
                )) %>%

        # Back transform log variables
        mutate(pred_wue = exp(pred_wue_log),
                pred_n_area = exp(pred_n_area_log),
                pred_pnue =  exp(pred_pnue_log),

            ) %>%

        # Remove log predictions and init height
        dplyr::select(-c(init_height, pred_wue_log, pred_n_area_log))
# Generate plot combinations

vars_q2_interaction <-

  crossing(

    # Get all numeric variables to plot (all y)
    as_tibble(t(combn(dplyr::select(data_pred_traits, where(is.numeric)) %>% names, 1))),

    # Select factor variables to plot
    x_axis_var = dplyr::select(data_pred_traits, nfixer) %>%  names,
    group_var = dplyr::select(data_pred_traits, treatment) %>%  names) %>%

    filter(V1 %in% c('pred_amax', 'pred_wue'))
vars_q2_interaction %>%
      # Gererate plots
      pmap( ~ boxplot_plot_pmap(data = data_pred_traits,
                                y = !!sym(..1), x = !!sym(..2),
                                fill = !!sym(..3)))
## [1] 24.08
## [1] 7.034
## [[1]]

## 
## [[2]]

vars_q2_treatment <-

  crossing(

    # Get all numeric variables to plot (all y)
    as_tibble(t(combn(dplyr::select(data_pred_traits, where(is.numeric)) %>% names, 1))),

    # Select factor variables to plot
    x_axis_var = dplyr::select(data_pred_traits, treatment) %>%  names,
    group_var = dplyr::select(data_pred_traits, treatment) %>%  names) %>%

    filter(V1 %in% c('pred_gs', 'pred_n_area', 'pred_pnue' ))
vars_q2_treatment %>%
      # Gererate plots
      pmap( ~ boxplot_plot_pmap(data = data_pred_traits,
                                y = !!sym(..1), x = !!sym(..2),
                                fill = !!sym(..3)))
## [1] 0.359
## [1] 3.399
## [1] 117
## [[1]]

## 
## [[2]]

## 
## [[3]]

Nodule weight

as_tibble(emmeans(models_q2$nodule_weight,
        pairwise ~ treatment,
        adjust = "tukey"
        )$contrast) %>%
        mutate(across(2:6, round, 6)) %>%
        kable()
contrast estimate SE df t.ratio p.value
ambientrain - ambientrain_nutrients -0.0631 0.0231 41 -2.7291 0.0443
ambientrain - ambientrain_water -0.0227 0.0129 41 -1.7523 0.3107
ambientrain - ambientrain_water_nutrients -0.0605 0.0198 41 -3.0512 0.0200
ambientrain_nutrients - ambientrain_water 0.0404 0.0239 41 1.6944 0.3397
ambientrain_nutrients - ambientrain_water_nutrients 0.0026 0.0265 41 0.0978 0.9997
ambientrain_water - ambientrain_water_nutrients -0.0378 0.0205 41 -1.8480 0.2663
# Treatment effects
emmeans_table_tidy(models_q2$nodule_weight,
                        formula = "treatment")
## [1] "Formula for pairwise comparisons: "
## pairwise ~ treatment
## <environment: 0x630b9a1089e8>

Nodule Count

as_tibble(emmeans(models_q2$nodule_count,
        pairwise ~ treatment,
        adjust = "tukey"
        )$contrast) %>%
        mutate(across(2:6, round, 6)) %>%
        kable()
contrast estimate SE df t.ratio p.value
ambientrain - ambientrain_nutrients -1.635 4.774 41 -0.3423 0.9860
ambientrain - ambientrain_water 18.157 4.774 41 3.8036 0.0025
ambientrain - ambientrain_water_nutrients -17.767 6.858 41 -2.5908 0.0611
ambientrain_nutrients - ambientrain_water 19.791 4.954 41 3.9951 0.0014
ambientrain_nutrients - ambientrain_water_nutrients -16.133 6.570 41 -2.4555 0.0825
ambientrain_water - ambientrain_water_nutrients -35.924 6.154 41 -5.8371 0.0000
# Treatment effects
emmeans_table_tidy(models_q2$nodule_count,
                        formula = "treatment")
## [1] "Formula for pairwise comparisons: "
## pairwise ~ treatment
## <environment: 0x630b9d2239d8>

Boxplots Nodules

data_for_nodule_predictions <- 
        data_nodules_cleaned %>%

        rownames_to_column("id") %>%

        # Remove unused variables
        dplyr::select(id, spcode, treatment, init_height)
string <- c("models_q2")

data_pred_nodules <-

        # Get models prediction
        gather_predictions(data_for_nodule_predictions,

                           # Return predictions
                            models_q2$nodule_count,
                            models_q2$nodule_weight) %>%    

        pivot_wider(names_from = model, values_from = pred) %>%
            rename_all(funs(

                # rename columns
                stringr::str_to_lower(.) %>%
                stringr::str_replace(., c(string),"pred_") %>%

                # Remove dollar sing
                gsub("\\$", "", .)
                )) %>% 
        # Remove log predictions and init height
        dplyr::select(-c(init_height))
# Generate plot combinations

vars_q2_nodules <-

  crossing(

    # Get all numeric variables to plot (all y)
    as_tibble(t(combn(dplyr::select(data_pred_nodules, where(is.numeric)) %>% names, 1))),

    # Select factor variables to plot
    x_axis_var = dplyr::select(data_pred_nodules, treatment) %>%  names,
    group_var = dplyr::select(data_pred_nodules, treatment) %>%  names) 
vars_q2_nodules %>%
      # Gererate plots
      pmap( ~ boxplot_plot_pmap(data = data_pred_nodules,
                                y = !!sym(..1), x = !!sym(..2),
                                fill = !!sym(..3)))
## [1] 146.5
## [1] 0.199
## [[1]]

## 
## [[2]]